Realistic Grand Canonical Monte Carlo Surface Simulation: Application to Ar(lll) 
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Most realistic, off-lattice surface simulations are done canonically — conserving particles. For some 
applications, however, such as studying the thermal behavior of rare gas solid surfaces, these con- 
stitute bad working conditions. Surface layer occupancies are believed to change with temperature, 
particularly at preroughening, and naturally call for a grand canonical approach, where particle 
number is controlled by a chemical potential. We report preliminary results of novel realistic grand 
canonical Monte Carlo simulations of the Lennard- Jones (LJ) fcc(lll) surface, believed to represent 
a quantitative model of e.g. Ar(lll). The results are successful and highly informative for temper- 
atures up to roughly 0.8 T m , where clear precursor signals of preroughening are found. At higher 
temperatures, convergence to equilibrium is hampered by large fluctuations. 

PACS numbers: 68.35.Rh, 68.35. Bs, 68.35.Md, 82.65.Dp 



1. Introduction 

Our quantitative understanding of thermal surface 
phenomena, including surface phase transitions, relies 
heavily on numerical simulations. Classical molecular dy- 
namics (MD) and Monte Carlo (MC) simulations have 
played and play a crucial role, describing the surface 
either via lattice hamiltonians, or via continuous (off- 
lattice) classical hamiltonians, or via continuous ab ini- 
tio hamiltonians. The continuous hamiltonians are of 
course much more realistic, but they also have serious 
limitations. One of these limitations, which we address 
here, has traditionally been connected with strict parti- 
cle conservation. It is desirable, in the study of many 
surface phenomena, to let the particle number fluctuate 
grand canonically. Lattice models permit that very nat- 
urally but continuous models do not. In fact, there is so 
far very little experience of grand canonical surface stud- 
ies conducted with continuous, off-lattice MD and MC 
simulation. In this paper we report on groundwork, plus 
some initial successes, in implementing a grand canonical 
Monte Carlo (GCMC) surface simulation based on a con- 
tinuous classical hamiltonian, in particular the Lennard- 
Jones fcc(lll) surface, chosen as a model for rare gas 
surface. 

Rare gas solid surfaces and films have provided an im- 
portant testing ground for several surface phase transi- 
tions for over two decades. Surface melting, roughening, 
and recently preroughening have been studied and char- 
acterized at the free solid- vapor interface. Layering tran- 
sitions of thin rare gas film on smooth substrates have 
given rise to a vast literature. The discovery of reen- 
trant layering of rare gases on substrates has led to a 



debate |l| over a possible explanation in terms of pre- 
roughening, as suggested by RSOS models ||, versus a 
melting-solidification interplay |3j as argued from canon- 
ical simulations with a continuous potential. In real life, 
the true rare gas solid surface begins to develop diffu- 
sion phenomena, which are totally absent in lattice RSOS 
models, above roughly (1/2) T m . However, the lattice 
models allow a much more thorough statistical mechan- 
ical understanding and classification of possible surface 
phases. In particular, they imply the possible existence 
of a preroughening (PR) transition, leading to a so called 
Disordered Ordered Flat (DOF) phase. At Tpr, among 
other things, the surface occupancy exhibits a transition 
from close to a full monolayer (T < Tpr) to close to a 
half monolayer (T > Tpr). This kind of transition can- 
not be easily simulated with a fixed particle number. In 
fact it has been shown elsewhere |^,|| that a fixed particle 
number will lead, at and above Tpr,, to a kind of phase 
separation into two neighboring DOF phases. In other 
words, that surface will not even remain macroscopically 
homogeneous. Classical off-lattice GCMC [|| has so far 
be succesfully applied to simulate different systems such 
as capillary condensation in nanopores j?J . However diffi- 
culties appear when one tries to simulate a denser system 
H and the classical GCMC technique requires modifica- 
tions. In the following we shall describe our own imple- 
mentation. 

2. Grand canonical Monte Carlo surface simula- 
tion: technical details. 

Our implementation of GCMC is based on well known 
techniques |1, with modifications to improve the sam- 
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pling efficiency in a system with a free surface. In partic- 
ular, taking inspiration from an idea previously used for 
bulk liquids [|| , we have restricted attempts to create or 
destroy particles to a region near the surface, where they 
are much more likely to be accepted. 

We have considered three different types of Monte 
Carlo moves: 

(a) small thermal displacements of individual atoms, 
with magnitudes Sr adjusted with T in order to obtain an 
optimal acceptance rate of 50%, and random direction; 

(b) large lateral displacements of surface atoms; 

(c) destruction of an existing particle, or creation of a 
new particle. Since deep in the bulk the acceptance rates 
will generally be exceedingly small, we restrict attempts 
to within a region of thickness d = na, where a is the 
spacing between two consecutive bulk lattice layers, and 
n is typically 6. This region is centered on the outer layer 
of the system simulated. 

Moves of type (b) were thought at first to be impor- 
tant for establishing surface diffusivity. It was later seen, 
however, that this is not really so for the LJ surface. 
Therefore, they were eventually omitted. Clearly, only 
moves of type (c) are grand canonical. Thermal moves 
(a) must be by far the most frequent, for equilibration. 
We find that a large number of such moves is necessary, 
in the temperature range investigated, between a grand 
canonical move (c) and the next. We have taken the 
probability of attempting a move (c) to be a c = ad = a 
(respectively for creation and destruction), and that of 
attempting a move (a) to be a m = 1 — a, with a very 
small, its actual value depending on the ratio between 
the number of atoms in the outer layers and the total 
number of atoms, but roughly of the order of 0.5 x 10~ 3 . 

The acceptance probability of creation and destruc- 
tion, p c and p ( i have then been finally normalized in the 
form H: 
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where N is the number of atoms present in the 
creation/destruction region, V its volume, A = 
[/i 2 /(27rmfcsT)] 1 / 2 is the thermal de Broglie wavelength, 
and SU — U a — Ub is the difference between the total 
energy U a after the trial move (creation or destruction) 
and the total energy Ub before it. In order to satisfy de- 
tailed balance N must be reajusted at each step of the 
simulation. 

We simulated in this way the free fcc(lll) solid-vapor 
interface of a rare gas, particularly Ar. The interatomic 
forces were described by a (12,6) Lennard-Jones poten- 
tial truncated at 2.5 er. Our system consisted of a 15-layer 
slab, with periodic boundary conditions along the x and 



y direction in the interface plane. Three bottom atomic 
layers were kept frozen, while 12 layers of 480 atoms each 
were free to evolve. The lateral box size was rigid, but 
readjusted at each temperature according to the ther- 
mal expansion coefficient of this potential, obtained from 
separate bulk simulations 9. We considered a grid of 
temperatures above 0.46e, as there is very little action 
below (the bulk melting temperature is Tb — 0.7e). For 
each temperature we first found the equilibrium value of 
the effective chemical potential fi a (T). This was done 
by trial and error, starting from an arbitrary value, and 
changing it until the particle number remained as sta- 
tionary as possible. With HoiT) so determined we made 
long equilibration runs looking for stabilization of both 
the total energy and the number of particles in the sys- 
tem. Generally half a million of Monte-Carlo (MC) steps 
(in the usual sense, i.e., one step attempting to move 
each particle on average once) were sufficient to reach 
reasonable equilibrium . Here we encountered our main 
problem, which is that this procedure is increasingly less 
convergent, and the resulting surface less stable, as T in- 
creases. In practice, we found it impossible to stabilize 
the surface against exceeedingly large fluctuations above 
T = 0.56e (~ 0.80 T m ). This constitutes an important 
drawback, because the interesting phase transitions of 
the surface are believed to lie just above this tempera- 
ture, i.e., preroughening at 0.83 T m and roughening at 
0.95 T m 0. On the other hand, a further reduction in 
the relative creation/destruction rate a, necessary to ob- 
tain grand-canonical equilibrium at these higher temper- 
atures, would quickly become too costly, as would the 
alternative possibility of creating/destructing atoms di- 
rectly in the vapor phase. Below 0.80 T m , where our equi- 
libration worked, typically 30 to 50 uncorrelated config- 
urations were subsequently generated from another half 
a million MC steps run, and finally analysed. 

3. Results and discussion 

The typical (xy)-averaged density profile of the simu- 
lated grand canonical Ar(lll) surface (really solid- vapor 
interface, but the vapor density is tiny) is shown in fig. 
|l|. Layers are clearly visible, and the outermost ones 
are labeled 2 ( "adatom layer" ) , 1 ( "surface layer" ) , and 
("subsurface layer"). We conventionally associate an 
atom with a given layer n, if it falls within the bin n, as 
indicated. Taking the ratio of the number of atoms in 
layer n, N n to the number Nb of atoms in a bulk layer, 
we obtain a conventional layer occupancy o„ = N u /Nb. 
Occupancies of the three outermost layers are shown in 
fig. H as a function of T. At low temperatures, the layers 
and 1 were almost full, with only a few percent of va- 
cancies, with corresponding few adatoms in layer 2. As 
T increased, layer 1 lost atoms very rapidly, reaching an 
occupancy o\ ~ 55% at T — 0.56e. 
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FIG. 1. (xj/)-averaged density profile of the simulated Ar 
fcc(lll) surface in grand canonical equilibrium at T — 0.56e, 
or T = 0.8 T m , plotted along the surface normal. The 
sub-monolayer coverage in layer 1 is the main new result, in 
comparison with previous canonical simulations. 



canonically o 2 ~ 100%, Oi ~ 90%, oo — 10%, against 
the present grand canonical values 02 — 94%, 01 ~ 55%, 
o ~ 8%. Grand canonically, the "total surface popula- 
tion" (o 2 + o\ + Oo — 1) appears to change continuously 
with temperature, going from nearly 100% at 0.65 T m 
down to 57% at 0.8 T m . This kind of gradual change was 
not predicted by the less realistic solid-on-solid models 




Extrapolating, these data strongly suggests that at 
T ~ 0.58e the surface spontaneously tends to half oc- 
cupancy, with about 8% adatoms and 6% subsurface va- 
cancies. This temperature is in excellent agreement with 
the experimentally observed temperature for the onset of 
reentrant layering, T = 69 K ~ 0.83 T m 0. The coinci- 
dence with the preroughening temperature obtained from 
canonical simulations is perfect. The half occupancy 
also agrees with theoretical predictions and with solid-on- 
solid simulations of the disordered flat (DOF) state, real- 
ized at and above preroughening [^,^2). Therefore, the 
present grand canonical realistic simulation data strongly 
confirm the occurrence of preroughening and of a DOF 
phase on the free surface of Ar(lll) at Tpr = 0.83 T m . 
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FIG. 2. Occupancies of the three outermost layers (frac- 
tion of particles contained in the bins indicated in fig. hi) as 
a function of temperature. Note the gradually decreasing oc- 
cupancy of layer 1, extrapolating to the preroughening value 
of 1/2 close to T = 0.58e. 



It is important at this point to compare the properties 
of this grand canonical Lennard- Jones (111) surface with 
those of exactly the same surface, simulated under identi- 
cal physical conditions except for the canonical, particle- 
conserving constraint. The results obtained in our ear- 
lier, very extensive MD simulations |^| were quite differ- 
ent. By starting with an integer number of layers, and at 
the very same T — 0.56e, for example, we had obtained 




FIG. 3. Top view (snapshot) of the three outer layers of 
simulated grand canonical Ar(lll) at 0.8 T m . This surface is 
nearly "disordered flat" (DOF). Atoms in layers 0, 1 and 2 , 
identified through their ^-coordinate as in fig. [I], are respec- 
tively black, grey and white. For convenience, four adjacent 
cells are shown. The nearly half occupancy of the surface layer 
(layer 1) is realized through large islands and large craters. 

Fig. ^| presents a top snapshot of the grand canonical 
surface at T = 0.56e, the closest we can get to Tpr. This 
is, we believe, the first available illustration of what a rare 
gas surface really looks like at preroughening. The main 
features to be noted are the large islands and craters in 
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the first (gray) layer. We did not try yet to examine the B 48, 14556 (1993). 

height-height correlations of this surface, which would [12] M. den Nijs and K. Rommelse, Phys. Rev. B 40, 4709 
have been very interesting, particularly to check whether (1989). 
they are large as expected by our predicted divergence 
at Tpr fl. In order to do that properly, we would need 
finite-size scaling, which is computationally too demand- 
ing. 

4. Conclusions 

In conclusion, we have reported preliminary results on 
the microscopic nature of a rare gas solid- vapor interface, 
as obtained by means of a novel grand canonical Monte 
Carlo simulation. In spite of difficulties encountered with 
stabilizing the surface at higher temperatures, we have 
succeeded in describing what appears to be a well equi- 
librated surface in an interesting temperature range. We 
have obtained a clear indication that preroughening is 
indeed taking place under conditions very close to those 
where reentrant layering was observed on Ar(lll), and 
where DOF phase separation was found by canonical MD 
B. We expect that, with the improvements which expe- 
rience should soon bring about, grand canonical surface 
simulations should become the standard approach in the 
near future. 
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